% This code estimates the entry cost distribution parameters for each
% market.

function [mu_f_all, sigma_f_all,...
    mean_cost_entry, std_cost_entry, entry_exitflag_vec]...
    = entry_cost_estimation(D_input, profit, N_st, N_f, estimate_mu)

% get data
major      = D_input.major;
n_insurers = D_input.n_insurers;
mkt        = D_input.mkt;
NN = length(major);

% for each market, compute the symmetric fringes' profit & # of fringes
profit_f_vec = zeros(N_st, 1);
n_f_vec = zeros(N_st,1);

for m=1:N_st
    % indices of fringes in mkt m
    ind_f_m = find(mkt==m & major==0); 
    
    % profit
    profit_f_m = profit(ind_f_m); 
    
    if length(unique(profit_f_m))==1
        profit_f_vec(m) = unique(profit_f_m);
    elseif length(unique(profit_f_m))>1
        if abs(max(unique(profit_f_m))-min(unique(profit_f_m)))<1e-10 % numerical approximation
            disp('numerical approximation: entry_cost_estimation.m')
            profit_f_vec(m) = mean(profit_f_m);
        else
            disp('ERROR1 in entry_cost_estimation.m') % diff too big
        end
    elseif length(unique(profit_f_m))<1 % each mkt has at least one fringe
        disp('ERROR2 in entry_cost_estimation.m')
    end

    % store the number of fringes in mkt m
    n_f_vec(m) = unique(n_insurers(ind_f_m));
end

% log normal entry cost parameters
% G(cost_entry_threshold) = n_f/N_f where N_f is the potential # of firms
if estimate_mu==0 % estimate sigma
    mu_f = 5;
else % estimate mu
    sigma_f = 0.6; 
end

if estimate_mu==0
    sigma_f_vec = zeros(N_st,1);
else
    mu_f_vec = zeros(N_st,1);
end

fval_vec = zeros(N_st,1);
entry_exitflag_vec = zeros(N_st,1);

options = optimoptions('fsolve','Display','none', ...
    'FunctionTolerance', 1e-15, 'OptimalityTolerance', 1e-15);

% for each mkt, compute the entry cost distribution
for m=1:N_st
    
    cost_entry_m = profit_f_vec(m); % cutoff entry cost is the per-fringe profit
    
    n_f_m = n_f_vec(m);
    
    if estimate_mu==0
        ff = @(sigma_f)solve_f_sigma(sigma_f, mu_f, cost_entry_m, n_f_m, N_f);
        x_ini = 0.5;
    else
        ff = @(mu_f)solve_f_mu(mu_f, sigma_f, cost_entry_m, n_f_m, N_f);
        x_ini = 1; 
    end

    [x, fval, exitflag] = fsolve(ff, x_ini, options);
    
    if estimate_mu==0
        sigma_f_vec(m) = x;
    else
        mu_f_vec(m) = x;
    end

    fval_vec(m) = fval;
    entry_exitflag_vec(m) = exitflag;
end

% display results
if estimate_mu==0
    [mean_cost_entry, var_cost_entry] = lognstat(mu_f, sigma_f_vec); % N.st x 1
else
    [mean_cost_entry, var_cost_entry] = lognstat(mu_f_vec, sigma_f);
end
std_cost_entry = sqrt(var_cost_entry); % N.st x 1

temp =  zeros(2,5);
temp(1,:) = [mean(mean_cost_entry(entry_exitflag_vec>0)) ...
    median(mean_cost_entry(entry_exitflag_vec>0)) ...
    min(mean_cost_entry(entry_exitflag_vec>0)) ...
    max(mean_cost_entry(entry_exitflag_vec>0)) ...
    std(mean_cost_entry(entry_exitflag_vec>0))];
temp(2,:) = [mean(std_cost_entry(entry_exitflag_vec>0))...
    median(std_cost_entry(entry_exitflag_vec>0))...
    min(std_cost_entry(entry_exitflag_vec>0))...
    max(std_cost_entry(entry_exitflag_vec>0))...
    std(std_cost_entry(entry_exitflag_vec>0))];
vartype = {'E(c)', 'Std(c)'};

disp(' ')
disp('****************')
disp('Market-specific entry cost')
disp('               Mean       Median       Min        Max         Std ')
for hh=1:2
    fprintf('%-8s %10.0f %10.0f %10.0f %10.0f %10.0f \n', ...
        vartype{hh}, temp(hh,:));
end

clear temp

if estimate_mu==1
    temp = [mean(mu_f_vec(entry_exitflag_vec>0)) ...
        median(mu_f_vec(entry_exitflag_vec>0)) ...
        min(mu_f_vec(entry_exitflag_vec>0)) ...
        max(mu_f_vec(entry_exitflag_vec>0)) ...
        std(mu_f_vec(entry_exitflag_vec>0))];
    disp(' ')
    disp('****************')
    disp('Market-specific mu')
    disp('      Mean       Median       Min        Max         Std ')
    fprintf('%10.2f %10.2f %10.2f %10.2f %10.2f \n', temp);
    disp(['sigma fixed to = ', num2str(sigma_f)])
else
    temp = [mean(sigma_f_vec(entry_exitflag_vec>0))...
        median(sigma_f_vec(entry_exitflag_vec>0))...
        min(sigma_f_vec(entry_exitflag_vec>0))...
        max(sigma_f_vec(entry_exitflag_vec>0))...
        std(sigma_f_vec(entry_exitflag_vec>0))];
    disp(' ')
    disp('****************')
    disp('Market-specific sigma')
    disp('      Mean       Median       Min        Max         Std ')
    fprintf('%10.2f %10.2f %10.2f %10.2f %10.2f \n', temp);
    disp(['mu fixed to = ', num2str(mu_f)])
end

% map estimated entry cost parameters to a vector of size NN x 1
mu_f_all = zeros(NN, 1);
sigma_f_all = zeros(NN,1);

for m=1:N_st

    ind_f_m = find(mkt==m & major==0); % indices of fringes in mkt m

    if estimate_mu==0
        mu_f_all(ind_f_m) = mu_f; % same for all mkts
        sigma_f_all(ind_f_m) = sigma_f_vec(m); 
    else
        mu_f_all(ind_f_m) = mu_f_vec(m);
        sigma_f_all(ind_f_m) = sigma_f; % same for all mkts
    end
end


